################################################################################
# Code for the analysis of the Replication
# Estimating Onsets of Binary Events in Panel Data
# Liam F. McGrath
#
# Purpose: This code analyses the results of the sensitivity analysis for
#          the determinants of the onset of civil war conducted in the paper
#          (and based upon Hegre and Sambanis 2006).
#          This script requires the output of the HS replication, created
#          by the script eba_hs_run.R 
#          If you do not wish to run the script yourself these outputs can be
#          found in the output subfolder
################################################################################

rm(list = ls())

library(ggplot2)
library(ggthemes)
library(plyr)
library(reshape2)
library(xtable)


# IMPORTANT:
# Define the path which all the code, data and results are located within
# From this changes in working directory are formed by pasting this
# string with a string for the relevant subfolder.
# If using this code be sure to change to this to your own directory
dir <- "/Users/liammcgrath/Documents/Academic/Research/transition/replication_materials/"

setwd(paste0(dir, "output"))

load("summaries.RData")

############
# TABLE OF COEFFICIENTS (Table 5 in the Appendix)
############



# Table with vars ranked by p-value for effect in models where ongoing years
# are set to zero
# Focus on those weighted by log-likelihood

rank_table <- summaries[, c("varname", "varconc", "beta_w_zero", "beta_w_na",
                            "beta_w_mark","pvalue_w_zero", "pvalue_w_na", 
                            "pvalue_w_mark")]

rank_table <- rank_table[order(abs(rank_table$pvalue_w_zero)),]

#  print(xtable(rank_table), include.rownames = FALSE, file = "hs06_rep_full.tex")
# 


##########
# Constructing the Graphs for the paper
##########

# To simplify, focus on cases where weighted p-value less than 0.05
# for at least one model

p_comp <- 0.05 # pvalue to categorise effects for comparison


# Compare NA to zero
s_na_s_zero <- summaries[summaries$pvalue_w_na <= p_comp & 
                           summaries$pvalue_w_zero <= p_comp, 1]


ns_na_s_zero <- summaries[summaries$pvalue_w_na > p_comp & 
                            summaries$pvalue_w_zero <= p_comp, 1]


s_na_ns_zero <- summaries[summaries$pvalue_w_na <= p_comp & 
                            summaries$pvalue_w_zero > p_comp, 1]


ns_na_ns_zero <- summaries[summaries$pvalue_w_na > p_comp & 
                             summaries$pvalue_w_zero > p_comp, 1]

#######
# Ropeladder graphs Included
#####

# Create a variable id ordered by p-value when ongoing set to zero
rank_table$varname_i <- factor(rank_table$varname, 
                               levels = rank_table[order(-rank_table$pvalue_w_zero), "varname"])

rank_long_p <- melt(rank_table[,c("varname_i","pvalue_w_zero", "pvalue_w_na",
                                  "pvalue_w_mark")], id=c("varname_i"))


##########################################
# FIGURE 6: Results of replication of HS #
##########################################

# First define a vector with more informative variable names
varlabels <- c("Neighbour at War", "Partially free polity", "Rough terrain",
               "GDP Growth", "Military personnel (in thousands)",
               "Years Since Last Regime Change (decay function)",
               "Oil Exports as a proportion of GDP", 
               "Middle East and North Africa Region Dummy",
               "Regulation of Participation",
               "Political Instability")

# Plot of p-values (right panel)
p_plot <- qplot(data = rank_long_p[rank_long_p$varname_i %in% c(ns_na_s_zero, s_na_ns_zero, s_na_s_zero),],
                x = varname_i, y = value, group = variable, shape = variable, 
                size = I(3.5)) + 
                coord_flip() + 
                geom_hline(yintercept = 0.05, linetype = "dashed") +
                ylab("Non-normal weighted p-value")  + theme_bw() + 
                scale_x_discrete(labels = NULL) +  # remove labels for sake of format of paper
                xlab(" ") +
                scale_shape_discrete(name="Model", 
                        breaks=c("pvalue_w_zero", "pvalue_w_na", "pvalue_w_mark"), 
                        labels=c("Zero", "Missing", "Markov"))
# ggsave("pvalue_rope.pdf", p_plot, width = 7, height = 7)

# Plot of relative change in Beta (left panel)
rank_b <- ddply(rank_table, .(varname_i), summarise,
                beta_zero = rep(0, length(varname_i)),
                beta_na = (beta_w_na - beta_w_zero) / beta_w_zero,
                beta_mark = (beta_w_mark - beta_w_zero) / beta_w_zero
                )

rank_long_b <- melt(rank_b, id=c("varname_i"))
b_plot <- qplot(data = rank_long_b[rank_long_b$varname_i %in% c(ns_na_s_zero, s_na_ns_zero, s_na_s_zero),],
                x = varname_i, y = value * 100, group = variable, 
                shape = variable, size = I(3.5)) + 
                coord_flip() + geom_hline(yintercept = 0, linetype = "dashed") +
                ylab(expression(paste("Weighted Mean ", beta, " relative to ", beta,
                        " when ongoing years set to zero (%)"))) +
                xlab("Variable") + theme_bw() + theme(legend.position="none") +
                scale_x_discrete(labels= varlabels) 
# ggsave("beta_rope.pdf", b_plot, width = 9, height = 7)
